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Abstract. A new approach to the kinetic simulation of plasmas in complex 
geometries, based on the Particle-in- Cell (PIC) simulation method, is explored. In 
the two dimensional (2d) electrostatic version of our method, called the Arbitrary 
Curvilinear Coordinate PIC (ACC-PIC) method, all essential PIC operations arc 
carried out in 2d on a uniform grid on the unit square logical domain, and mapped 
to a nonuniform boundary-fitted grid on the physical domain. As the resulting logical 
grid equations of motion are not separable, we have developed an extension of the 
semi-implicit Modified Leapfrog (ML) integration technique to preserve the symplectic 
nature of the logical grid particle mover. A generalized, curvilinear coordinate 
formulation of Poisson's equations to solve for the electrostatic fields on the uniform 
logical grid is also developed. By our formulation, we compute the plasma charge 
density on the logical grid based on the particles' positions on the logical domain. 
That is, the plasma particles are weighted to the uniform logical grid and the self- 
consistent mean electrostatic fields obtained from the solution of the logical grid 
Poisson equation are interpolated to the particle positions on the logical grid. This 
process eliminates the complexity associated with the weighting and interpolation 
processes on the nonuniform physical grid and allows us to run the PIC method on 
arbitrary boundary-fitted meshes. 



PACS numbers: 52.20.Dq,52.25.Dg,52.35.Fp,52.65.Rr 
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1. INTRODUCTION 

The standard Particle-in-Cell (PIC) method utilizes macroparticles to follow the phase 
space evolution of weakly coupled (collisionless) plasmas. In this regime, PIC can be 
thought of as an approximate method of integration for the Vlasov equation pQ: 

^7 + v ■ V g f s + —(E + vxB)- Vefs = 0, (1) 
at m s 

where f s (x, v, t) is the time-dependent six- dimensional phase space distribution function 
of the plasma particles of species s. For simplicity, in this work we limit ourselves to 
the electrostatic case with no applied background magnetic field. Applying the method 
of characteristics [2] to ((H), it can therefore be shown that the PIC macroparticle orbits 
evolve self-consistently in phase-space according to 

? = * (2) 

$ = _^ALV$(f). 1 ' 

In 02]), an overdot denotes a time derivative and the macroparticle charge and mass 
satisfy ^m/^m = ?sA"s- The mean- field electrostatic potential is obtained on a 
computational mesh via Poisson's equation 

n m 

V 2 $(f) = -4vrp(f) = -4tt QmS(x - x,) . (3) 

i=l 

In our notation, Nm is the number of macroparticles (assumed to be of the same 
species with equal charge for simplicity) and Xi are the particle positions. p(x) is the 
charge density obtained at discrete locations on the mesh from the macroparticles at 
Xi using interpolation functions, S(x — xi), typically chosen to be _B-splines j3]. These 
interpolation functions effectively give the macroparticles a finite width based upon the 
grid spacing and lead to cutoff Coulombic interactions between macroparticles [I]. 

PIC codes are generally designed using rectangular meshes in Cartesian geometry, 
but have been extended to cylindrical and spherical coordinates. However, extension to 
arbitrary grids has proven much more difficult. Jones [5] was among the first to develop a 
curvilinear-coordinate PIC method capable of operating on boundary-conforming grids 
tailored to accelerator and pulsed-power applications. Other codes were soon developed 
in an effort to model ion diodes [HI [7] and microwave devices [8] more accurately. 
These early methods involved generating a nonuniform initial grid based upon the 
physical boundaries of the system and running the PIC components on this physical 
grid. There are many benefits to this type of system, such as having smoothly curved 
boundaries in contrast to the "stair-stepped" boundaries inherent to the rectangular- 
grid PIC approach, which occur when a part of the boundary is not aligned with either 
coordinate surface. Furthermore, higher grid density can be placed in areas of interest 
within the system either statically or by allowing the grid to adapt dynamically [H [10] 
to the problem by following a pre-specified control function. Implemented wisely, such 
techniques should allow complex geometries to be simulated at a fraction of the cost 
associated with using a uniform grid code, in which the entire mesh must have the 
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resolution required to resolve the smallest physical features of the system, and in which 
the grid does not conform to the boundary. 

Arbitrary grid methods are not without their problems. Nonuniform grid cells 
make it computationally expensive to locate macroparticles on the grid for the charge 
accumulation and field interpolation steps of the PIC method, as well as for enforcing the 
particle boundary conditions. In a standard PIC code with a uniform structured grid, 
particle positions on the grid are easily determined. On a nonuniform grid, this location 
must be done iteratively [XT], which increases the computational cost of the method. 
Furthermore, interpolation on a nonuniform grid is complicated by the variations of 
the grid, which change the shape of the interpolation functions, often in a non-trivial 
way [12]. Finally, as the ratio of the largest to the smallest cell size increases and the 
number of particles per cell in the smallest cells becomes small, the amount of noise near 
the small cells also increases (assuming the charge density is roughly constant). Thus, 
complex and often time-consuming gridding strategies [13] and/or particle splitting and 
merging algorithms [2] must be implemented to keep the noise within the system to 
a minimal threshold value while the structures of interest within the system are still 
resolved. 

We consider the problems associated with these and other existing adaptive grid 
PIC approaches to be serious. As such, we have developed a new nonuniform grid 
PIC method, which incorporates some of the best features of several existing methods 
with a new idea for the implementation of the PIC method. Our goal is to design an 
arbitrary, curvilinear-coordinate PIC (ACC-PIC) code capable of operating efficiently 
and accurately on an arbitrary (but structured) moving mesh for a boundary of arbitrary 
geometry. We construct a logical (or computational) grid on the unit square and map it 
to the physical domain, as illustrated in Figured] We implement the main components of 
the PIC method-the charge accumulation, particle push, field solve, and interpolation- 
on this logical domain. This approach deals with all the problems listed in the previous 
paragraph except the last. The issue of having some cells with much fewer macroparticles 
than others is still a problem requiring attention, for example particle splitting and 
merging. 

In this paper, we present results on the development of these methods into a 2D, 
electrostatic PIC code. By implementing the PIC components on the logical grid, we 
show that we can eliminate the particle location and interpolation problems that plagued 
the earlier boundary conforming non-uniform grid methods [SI El [7J [8] . Particle locations 
are easily found on the logical grid using the same techniques as standard uniform grid 
PIC codes. Since the charge accumulation and field interpolation are both done on 
the logical grid, we are again able to use the efficient algorithms that are utilized in a 
uniform grid code. This eliminates the need for calculating non-standard particle shapes 
on the physical grid [12]. Furthermore, we develop a Hamiltonian-based, semi-implicit 
second-order accurate symplectic logical grid mover and apply it to the time advance of 
the particles that includes the effects of inertial forces. This mover has the important 
property of requiring only a single field solve per timestep. Since we are moving particles 
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Figure 1. Schematic for mapping from the logical to the physical grid for an annulus. 
The numbers along the boundaries of the square logical domain are used to indicate 
which edge of the square maps to which boundary segment on the physical domain. 



on a square grid, particle boundary conditions, which may be difficult to implement on 
a nonuniform physical grid, are fairly simple with our approach. Finally, our method 
allows us to solve the electrostatic field equations on a simple square mesh in the logical 
domain rather than in the complex physical domain. 

The remainder of this paper is structured as follows. In Section [21 we quickly review 
the Winslow method of generating a grid that conforms to a curved boundary, with finer 
grid near parts of the boundary with more curvature. We derive the particle equations 
of motion on the logical grid in Section [31 and discuss the semi-implicit modified 
leapfrog integration technique utilized for the integration of these equations in Section HJ 
Section [5] describes the solution of the electrostatic field equations on the logical grid, as 
well as details of the charge weighting/field interpolation schemes utilized. We perform 
several standard tests for the verification of our 2D, electrostatic, nonuniform grid PIC 
code in Section [7J Finally, we present our conclusions and suggestions for future work 
on the method in Section [SI A general overview of our differential geometry notation is 
given in Appendix A and a derivation of the Poisson equation in logical space is given 



in |Appendix B 



2. Logical to Physical Grid Mapping 

2.1. Grid Generation Technique 

For completeness, we include here a discussion of Winslow's Laplace method [15], a 
simple grid generation technique which illustrates the logical to physical grid mapping 
inherent to this work. In Winslow's method, a set of uncoupled Laplace equations are 
solved on a uniform square 2D logical domain, £,77 G [0 : 1]. With Winslow's method, 
the finest gridding is concentrated in the regions of highest boundary curvature, e.g. 
around an object at the center of the domain. (See Figure [H in which the surface 
labeled "1" might represent a dust grain, "3" represents an artificial boundary far away, 
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and symmetry is imposed on surfaces "2" and "4." As discussed in Section 5, this figure 
can represent either azimuthal or axial symmetry.) While this gridding method is rather 
primitive in that we have no control over where the grid is concentrated away from the 
boundary, we have chosen this method for its simplicity and suitability for illustrating 
the feasibility of our ACC-PIC method. Winslow's method provides a simple way of 
generating an initial boundary conforming grid, to be adapted to the solution at later 
timesteps, while still retaining conformation to the boundary using more sophisticated 
techniques. 

The Laplace equation can be written in any coordinate system, but for simplicity 
we have chosen to implement Winslow's method in the Cartesian coordinate system in 
physical space. In the physical domain, Laplace's equation takes the form 

8 8f a 

v ^ = a^ = °< <^ = v-,». <*> 

In 01]), the logical variables £ a are the dependent variables and the physical variables x 13 
are the independent variables. Here and in the rest of this paper we assume summation 
over repeated indices. However, since we would prefer to solve the Laplace equation 
on the uniform logical space (£, rj) rather than directly gridding (j4j) on the physical 
domain, solving for £(x,y),r)(x,y) and inverting, we must transform the set of PDE's 
such that x 13 are the dependent variables and are the independent variables. This 
transformation has been outlined by Liseikin [16] ; a detailed derivation can be found in 
|17j . After much manipulation, we write the final system of equations to be solved as 

x o x 1 (9^ x r\ 

& V o d 2 y , d i y n \° I 

g22g^-2gi2Q^ + giig$ = 0, 

where 

-* dx^ dx^ 

9^) = d^Qf' n,u,j = l,---,n (6) 

is the covariant metric tensor (see Appendix A). While ©(a) and (b) appear to be 
the same equation for x(£,r)) and for y(£,rj), they actually possess opposite boundary 
conditions along each segment of the grid boundary: For any segment of the boundary, 
we have Neumann boundary conditions for £ and Dirichlet boundary conditions for rj 
or vice-versa, both expressed in terms of rj) and y(£,r)). Furthermore, from (jj]) and 
these boundary conditions, we know that £(x,y) and rj(x,y) are conjugate harmonic 
functions. Thus, Winslow's method forces the grid lines to be orthogonal (i.e. gi2 = 0) 
in the physical space. 

Since the covariant terms g a p depend on ||^ and the solution we seek is 
x(£, r)),y(£, rj), the system of equations is nonlinear. (P§) are therefore discretized using a 
second-order accurate finite-difference method and, because of their non-linear form, are 
iteratively solved using an inexact Newton-Krylov solver with the Generalized Minimal 
Residual Method (GMRES) PI US]. 

Winslow grids can be specified for a wide array of physical domains by specifying 
boundary conditions on boundary segments, as illustrated in Figure [1] The grids 
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(c) (d) 

Figure 2. Grids generated using Winslow's method with a circular object along the 
inner physical grid boundary, symmetry along the bottom segments (and on the left in 
(d)), and an artificial boundary far away at top. Here we have used a 32 x 32 grid in 
order to show the mapping more clearly. The appearance of nonorthogonal grid lines 
is due to the straight lines used by the plotting tool. 



presented in Figure [2] are of some of the physical grids with circular objects along 
the inner boundary which we have generated using Winslow's method. As we discuss 
in Section 5, this circular segment can represent a cylindrical or spherical surface. We 
have also generated grids with elliptically shaped objects, and are capable of generating 
objects of fairly general shape and size on any boundary segment |17] . 



3. DEVELOPMENT OF LOGICAL GRID EQUATIONS OF MOTION 

The most obvious method for implementing our mover in logical space would seem to 
be by simply converting the Newton-Lorentz equations of motion from the physical to 
the logical grid. Specializing to Id for simplicity and defining u = £ = v/J, we recast 
(ED as 

£ = /(«) = u m 
* = 9iM = 1 J 

where J = ^| in Id and the prime symbol represents a derivative with respect to 
the logical coordinate £. For compactness we have dropped the subscripts designating 
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macroparticle quantities here and in the rest of this paper. Taking the formal divergence 
of ([7]), we have 



df 



dg_ 

du 



2J'u , , . 

+ 0. (8) 



C J 

The Newton-Lorentz equations of motion are therefore not divergence-free in these 
variables. As such, a time integration of these equations of motion with a standard 
"naive" leapfrog (LF) integrator, 

i' = £ + Af/(u) 
v! = u + Atg(£',u), 

will not preserve phase space area, and the particle orbits will typically spiral inward or 
outward with time [TTj . 

3.1. Hamiltonian Approach to Logical Grid Equations of Motion 

Clearly, since the system of equations in (JJJ) are not in a canonical Hamiltonian 
formulation, transformation of the Newton-Lorentz equations of motion into logical 
coordinates leads to a lack of (£, u) phase space area conservation. We therefore 
construct logical grid particle equations of motion based on Hamilton's equations by 
using a canonical transformation (x,p) — > (£, -P), where p a = mv a is the physical space 
momentum and P a is the logical space momentum. This transformation is constructed 
via an F 2 (x,P,t) generating function [20], thereby assuring Hamiltonian equations of 
motion on the logical grid. Under the canonical transformation we have 



P a = £;F 2 (x,P,t) 
r = Q§^F 2 (x,P,t). 



T^/T/T,' ( 10 ) 



We have chosen the F 2 generating function such that x and P are considered independent 
variables. The transformed Hamiltonian K is given by 

* = <»> 

A time-adaptive grid is beyond the scope of this paper; we therefore drop the second 
term of ( FlTT) . Specializing to the contact transformation £ = £(x), we write the F 2 (x, P) 
generating function as 

F 2 (f,P)=^(f)P /3 , (12) 

such that by ([lOD, we have p a = k^P 13 and £ a = £ a (f). Here k^ a is the inverse Jacobi 
matrix as defined in Appendix A. 

In the absence of a static background magnetic field B, the physical grid 
Hamiltonian is given in general form by 

^ = ES + 9$(f) ' (13) 

i 

where the sum is over three rectangular components in physical space. Notice that 
( TT5|) is fully separable, H = T{p) + V(x). In this work, we consider two distinct 2D 
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cases, with azimuthal (z) symmetry and axial ((f)) symmetry. In the former, the grid 
segment labeled "1" along the curved grid segment in Figure [1] represents, for example, 
a cylindrical probe. For the axisymmetric case, this segment would represent a spherical 
probe. We stress here that both cases lead to separable Hamiltonians in the physical 
domain. 

f|T3|) can then be used in conjunction with ( [IT]) to construct the logical grid 
Hamiltonian: 

K=-^(g» f P f, F')+V&, (14) 



p„ = _j_a^p /3p7 _^ = ^p ) . ( 15 ) 



where g 131 is the contravariant metric tensor as defined in (1A.8j) . Again, since 
gP J = g^iCr]), we must represent g^ 1 in terms of the covariant metric tensor as in 
( lA.lOj) . Note that on the logical grid the transformed Hamiltonian can be written 
K = T(£, P) + V(£), meaning that we have transformed the separable physical grid 
Hamitonian to an equivalent, but non-separable system. Applying Hamilton's equations, 
£ a = Jp^ and P Q = — J|| to (JHj) gives the logical grid equations of motion: 

We note here that the term quadratic in P represents the inertial force in these 
coordinates, and that both £ M and P M are functions of £ a and P a . Since these equations 
are obtained from a Hamiltonian in canonical variables, the divergence for each degree 
of freedom + |^ (no sum) is zero. 

4. PARTICLE PUSH: MODIFIED LEAPFROG INTEGRATOR 

Whereas the physical space Hamiltonian leads to separable equations of motion which 
can be integrated by standard LF integration techniques, in the logical space this is no 
longer true. Since ( fT5|) is not a separable system of equations, integration with the LF 
method will not conserve phase space area. 

As such, we have chosen to implement an extension of the semi-implicit modified 
leapfrog (ML) integrator originally developed by Finn and Chacon [2T] for integrating 
2D solenoidal flows in fluid dynamics and magnetic field lines in MHD codes. (Because 
the dependence of the Hamiltonian on the canonical momentum is given analytically, 
the exactly divergence free interpolations used in [21] are not necessary.) Rewriting (|T5l) 

as P = W and £ = U, respectively, and denoting explicit (implicit) updates with a 
superscript e (i), the ML integrator can be written as M^t = P^t ° £Xt> wnere 

S=f + Aftf(S,P) pe . / £' = 6 

P 1 = P ' 'I P' = Pi + At V(it, A) 1 1 



Combining, we have 

i' = i+Atu{i',p) 

P' = P + AtV(£',P). 



(17) 
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The map is implicit and must be done by means of Newton or Picard iterations. 
The map P^t * s explicit, and can therefore be applied directly. It can be shown that, 
unlike the standard LF scheme commonly utilized in PIC codes, this non-time-centered 
formulation of the ML scheme results in only first-order accuracy in At. To achieve 
second-order accuracy in time, we simply symmetrize the ML scheme by composing 
^At/2°-^At/2 ^At/2 ^Ai/2- As s h° wn in [21], the implicit-followed-by-explicit ordering in 
each pair of mappings retains the area preserving nature of the integrator in a 2D phase 
space (one degree of freedom.) This integrator is seen to be symplectic for arbitrary 
degrees of freedom because it can be derived from a generating function [20] 

F 3 (i' , P) = -CP a + AtK(C', P), (18) 

with 



p'° = -^F 3 (i',P) 

r = -g^F 3 (£',P). 



(19) 



The second half of the ML update (-P^t/2 °£kt/2) * s described by an F 2 (£, P') generating 
function. The alternation of the steps in £ and P gives second-order accuracy in At. 
The logical flow of the ML integrator as implemented on ( II 5p is 

CAt/2 ~> P'kt/2 ~* P'kt/2 Ckt/2- (20) 

We note here that the charge density is accumulated on the grid and the mean field 
solve is performed after the implicit position update step of the symmetrized ML mover. 
While this mover requires us to pass through the particle array twice per timestep, it 
allows us to accumulate the charge density and solve for the fields only once per timestep. 

With the exception of the case in Sec. 7.4, the boundary conditions we use are 
periodic, and we therefore do not need to describe particle reflection in the logical 
space. Reflecting conditions on particles can be implemented with some care in logical 
coordinates. 



5. FIELD SOLVER: GENERALIZED POISSON EQUATION 

In order to solve the electrostatic field equation on the logical grid, we first write the 
Poisson equation on the physical grid 

1 19 9$ 

jV-/V4 = 7 ^./^ = -4^, (2!) 

where p x is the physical charge density. Here, / is a geometry factor. For azimuthal 
symmetry, / = 1 and the Poisson equation takes the usual form 

V 2 $ = -4vrp x . (22a) 

For axisymmetry, we have f = r and the Poisson equation takes the form 

V • (rV$) = -Anrp x . (22b) 
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The generalized curvilinear coordinate formulation of Poisson's equation on the logical 
grid takes the form 

1 9 = "4^, (23) 



as derived in Appendix B. The logical density p^ is equal to Jp x , leading to 

W^Hl = -47T//. (24) 



In this form, we accumulate the logical density p^ and solve the Poisson equation by 
conservative differencing, as described in the next section. This logical grid solver 
removes the necessity of writing and maintaining multiple complex-geometry Poisson 
solvers for the various coordinate systems we wish to model. 

6. NUMERICAL IMPLEMENTATION 

For our 2D code, we have a wide range of grid choices on which to perform validation 
tests for our method in this and the next section. For the sake of brevity, we outline 
only two here, both of which have been designed such that analytical expressions for the 
metric tensors are easily obtained. In addition to the annulus grid generated numerically 
by the methods described in Section [2J we define a doubly periodic, nonuniform, 
orthogonal grid using: 

X ^min H~ (x max 3'min)(^ £g sin 27I"£) , . 

V = Vmin + (Umax ~ Vmin) (V + e g ^ n 2lT V) ■ 

Here £ mm , x max , ?/ mm , and y max are constants to scale boundaries of the physical grid to 
form a rectangle of arbitrary size. Here, e g is the nonuniformity parameter which controls 
the amount of nonuniformity of the grids. A doubly-periodic grid has been chosen for 
the implementation of the periodic field and particle boundary conditions utilized for 
the tests in the Section [0 We have also designed a doubly-periodic, nonorthogonal grid 
given by 

X = Xmin + Omax ~ ^min) (f + tg SU1 2tt£ SU1 2tT7]) 
y = ymin + Omax ~ J/min) (V + e g ^ n 2 7t£ sin 2^7]) 

to test the effects of non-zero cross terms (g 12 ) in our code, both in the field solver 
and the particle mover. For simplicity, we have constrained the grid nonuniformity 
parameter e g in (125aj) and (125bj) to be the same in each dimension. Notice that both 
(I25al) and (125bl) require that e g < j- so that the grid does not fold. 

6.1. Discretization of Particle Equations 

In two spatial dimensions, the logical grid particle equations of motion ( ([TBI) ) can be 
written as 

f = ±(g^p + g^p) 

ml ? V> (26a) 

V = i(g l2 P, + g 22 P v ), 
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and 

P -i(p2da 11 I op P a.9 12 I P 2dg 22 \ av(|) 

'/ 2m \ c or] 5 '/ or] V or] I or) 



(26b) 



Note that we have written (1261) in terms of the contravariant metric tensor, g^ v = g^ u (x), 
which is easily obtained as the inverse of the covariant metric tensor, g^ u {^ ) as in flA.lOj) . 

The third momentum component is generated at t — and held as a constant as 
the simulation progresses. This term contributes to the simulation through its inclusion 
in the effective potential term, V(£ ). In azimuthal symmetry, (xi,x 2 ) = (x,y) and the 
effective potential is V az i(£ ) = ^± + (7$(£ ), such that the ignorable-direction momentum 

does not contribute to the momentum update equations. However, for an axisymmetric 

p 2 

problem, (x\, 22) = {r, z) and the effective potential is V r ax i({ ) = 2mr ^ 2 +9^(0: suc ^ 
that its derivative is 

Of mr° ^ s 

(27) 



dv(j) niy% 

dr\ mr 3 



where {x^Xz) = (r, z) implies ju = dr/dC, and ju = dr/drj. Thus, in the axisymmetric 
case we must also interpolate the ju and j'12 components of the Jacobi matrix and 
the r-coordinate of the particle's physical space position to the particle position on the 
logical grid. 

6.2. Conservative discretization of Poisson equation 

In two dimensions the logical grid Poisson equation, (|24j) . takes the form 

d ( „a$ 19 <9$\ d ( 19 <9$ 99 <9$\ , , 

* \ D m + D \) + 1* [ D "* ikj] ~ (28) 

where we have defined D^ u = fjg^, again written in terms of the contravariant metric 
tensor g^ v {x). (12"H|) is comprised of a set of co-directed derivatives proportional to the 
diagnonal metric tensor components (D u and D 22 ) and a set of cross-directed derivatives 
for the off-diagonal metric tensor components (-D 12 ), i.e. V 2 $ = (V 2 $) co + (V 2 $) cross , 
the latter of which are nonzero for non-orthogonal coordinates (g 12 ^ 0). We have 
defined the coordinates £,77 at vertices thus $ and p are naturally defined at cell 
centers 2+1/2, j+1/2. The co- and cross-directed terms are discretized using appropriate 
centered-difference schemes, leading to the final discretized form of flUJ) : 

-4nfp r 



22 / /f> . „ _fls . . I _ r>22 



D 11 1 (i 



1 



4 A£ Ar/ 



Art 2 

D} 2 ( $, , 3 , , 3 - *. , 1 , , 1 ) - D} 2 $, 
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Note that, for a uniform grid with x, y G [0 : 1], the D 11 and .D 22 terms become 
unity (assuming the geometry factor / = 1) and D 12 = 0, thus f l29|) reduces to the 
standard 5-point discretization commonly used in Cartesian PIC codes. 

We note that this same discretization of (I28|) can also be obtained by the 
minimization of the variational principle [22] 



W 



IV$| 2 



dV, (30) 



discretized on the logical grid. The variational principle approach guarantees that, 
for properly applied boundary conditions, the matrix formed by the application of the 
discrete form of the V 2 operator is symmetric (for appropriate boundary conditions) and 
negative definite. This property is important since it permits the use a fast conjugate 
gradient (CG) solver. Symmetry is important because CG requires a symmetric, positive 
(or negative) definite matrix to converge to the correct solution, but for a non-symmetric 
matrix it may converge to the wrong answer. 

6.2.1. Validation of Poisson Solver using the Method of Manufactured Solutions The 
Method of Manufactured Solutions (MMS) [23] technique is utilized to validate our 
Poisson solver. MMS is a simple, yet powerful tool to construct solutions to PDE 
problems. By obtaining analytic solutions to problems, we are able to check that the 
error in the numerical solution converges to the analytic solution with second-order 
accuracy in grid spacing. MMS is based upon choosing a potential $ which satisfies a 
given set of boundary conditions a priori, then taking the required derivatives to find 
the source term, i.e. the density. To validate our solution of Poisson's equation, we 
simply choose a potential that satisfies a chosen set of boundary conditions, calculate 
the charge density, and use that density in our solver. Upon convergence of the solver, 
the numerical solution and the exact solution are compared, and the error between the 
two is calculated using the L 2 -norm: 



'V^f ((f) nnm — (f,MMS\ 2 

J ^ l=1 y i i — '-. (31) 

To ensure the accuracy of our method, we have set up MMS tests appropriate for both 
the analytically given grids of (25) as well as for the case in which an annular grid is 
numerically generated using the techniques of Section [2j 

For a unit square physical domain, we have chosen an MMS potential given by 

$mms = sin(27rx) sin(27ry), (32) 

such that, assuming Cartesian geometry (/ = 1), the MMS charge density is 

Pmms = 27r sin(2vrx) sin(27ry) . (33) 

For this particular choice of $mms, we test our Poisson solver with periodic field 
boundary conditions on all boundary segments. Table CD shows the L 2 -norm of the 
error in ($ — $mms) as given by ( 13TT) for the orthogonal, nonuniform grid ( (125a)) ) for 
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4.20043E-5 


1.47031E-4 


5.47060E-4 


128 


2.52982E-5 


1.19786E-5 


1.04332E-5 


3.65026E-5 


1.35787E-4 


256 


6.29958E-6 


2.98229E-6 


2.59901E-6 


9.09201E-6 


3.38199E-5 



Table 1. L 2 - n orm error between computational and analytic MMS potentials for 
different grid resolutions and grid nonuniformities (e ff ) for the 2D orthogonal grid 
given by (|25a[) . For these tests, we have chosen N% — N„. 

various levels of nonuniformity, e g . The results scale with second-order accuracy in grid 
spacing as expected, even for the e g = 0.15 case, in which the ratio of the area of the 

largest to the smallest cells is J mKX / = ("jrf^T") ~ 1140! 




e 

g 



Figure 3. Graph of the maximum skewness parameter S as defined in (|35|) as a 
function of e g for the grid given by (|25bj) . 





e g = 0.025 


€g = 0.05 


e g = 0.075 


€g = 0.15 


16 


2.08019E-3 


3.27156E-3 


5.57519E-3 


8.99093E-3 


32 


5.07340E-4 


8.20596E-4 


1.44430E-3 


2.34622E-3 


64 


1.25103E-4 


2.03916E-4 


3.62134E-4 


5.90944E-4 


128 


3.10476E-5 


5.07074E-5 


9.02599E-5 


1.47505E-4 


256 


7.73260E-6 


1.26353E-5 


2.25042E-5 


3.67914E-5 



Table 2. L2-norm error between computational and analytic MMS potentials for 
different grid resolutions and grid nonuniformities (e g ) for the 2D nonorthogonal grid 
given by ()25b|) . For these tests, we have again chosen = N v . 



Turning our attention to the nonorthogonal grid ( I25b[) . we must now also consider 
the amount of "skewness" of the grid as a major factor in the difficulty associated 
with solving the curvilinear Poisson equation (128j) . In an effort to characterize the 
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amount of nonorthogonality, or "skewness," of our grids, we can use the Cauchy- 
Schwartz inequality [21] , which states that any two points p and q in n-space must satisfy 
p-q< \p\\q\- Writing the contravariant metric tensor in the form g a/3 (x) = V£ Q • V£^, 
the inequality can then be rewritten as 

V£- Vt/ < |V£||Vt/|. (34) 

Squaring both sides and rearranging allows us to define the local grid skewness factor 

(W-Vr?) 2 (a 12 ) 2 

where < S < 1. In the limit of S(£,rj) — > 0, the grid is orthogonal at (£, rj), whereas 
for S(£,r)) — > 1, the grid becomes singular, i.e. the grid cells become elongated and can 
fold. Furthermore, we can define the maximum grid skewness, 

S max = max [Sfe, rjj)], (36) 

as a single parameter to characterize the maximum nonorthogonality of the generated 
grid. Figure [3] shows S ma , x as a function of e g for the grid given by ( I25bl) . For e g > 0.125, 
S'max is very nearly unity, and thus the grid is almost singular at some point on the grid. 
As noted above, this grid folds at e g = ^- ~ 0.16. As S mgx — > 1, the stiffness of the 
Laplacian matrix becomes quite large and the field solver will not easily converge, if at 
all. 

The MMS tests performed on the nonorthogonal grid used the same MMS potential 
as was used for the orthogonal square grid (and again / = 1). Table [2] shows that the L 2 - 
norm of the error scales with second-order accuracy in A£ = At/, as expected. We note 
here that the ratio of maximum cell area to minimum for the grid defined in f l25b[) scales 
according to = |^^ g (~ 34 for e g = 0.15), meaning that the ratio of the largest 
cell area to the smallest is much smaller than that of the orthogonal, nonuniform grid 
case, but the skewness of the grid adds to the challenge and leads to errors comparable 
to those in Table 1. 

We have also designed an MMS test to validate our Poisson solver for the annular 
grid generated in Section |2j We simulate half the annulus and apply symmetry 
conditions at the bottom. We apply either Dirichlet or Neumann boundary conditions 
along the inner and outer boundary segments, and symmetry requires homogenous 
Neumann boundary conditions along 9 = 0, ir. If the boundary conditions consist 
of only Neumann and periodic segments and no Dirichlet segments, the range of the 
Laplacian operator consists of densities that have total charge exactly equal to zero. 
Furthermore, this corresponds to a null space in the Laplacian operator $ — > $ + const. 
We have found that if we assure that the total charge is equal to zero, then the conjugate 
gradient algorithm converges quickly. The additional constant potential does not affect 
the electric field. We have therefore constructed the following potential: 

$ (r, 6) = 1 - r 3 + (r - n) (r 2 - r) cos 9, (37) 
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R = 5 


R = 10 


/? = 20 


16 


4.86684E-3 


8.78720E-3 


1.35298E-2 


32 


1.18612E-3 


2.15682E-3 


3.35229E-3 


64 


2.92357E-4 


5.32567E-4 


8.29718E-4 


128 


7.25471E-5 


1.32213E-4 


2.06105E-4 


256 


1.80677E-5 


3.29310E-5 


5.13433E-5 



(b) Spherical Laplacian 





R = 5 


= 10 


/? = 20 


16 


3.84324E-3 


6.41377E-3 


9.23903E-003 


32 


9.43887E-4 


1.60080E-3 


2.35531E-003 


64 


2.33099E-4 


3.96923E-4 


5.87121E-004 


128 


5.78702E-5 


9.86410E-5 


1.46102E-004 


256 


1.44142E-5 


2.45755E-5 


3.64121E-005 



Table 3. L2-norm error between computational and analytic MMS potentials for 
different grid resolutions and ratios R = ra/ri for the annular grids generated in 
Section [2] for both Cartesian (a) and cylindrical (b) solvers. For these tests, we have 
again chosen N% = N v 



where we have used polar coordinates to express (137|) in the physical space, as it more 
naturally aligns with the annular grid case than the Cartesian notation we have been 
using until this point. At this point we should mention 

As our Poisson solver is designed to solve general geometries in the physical space, 
we have applied the MMS test to both the azimuthally symmetric and axisymmetric 
cases (/ = 1 and f = r, respectively.) The physical Laplacian in this geometry with 
/ = 1, a cylindrical annulus, is given by (I22aj) . leading to the physical charge density 



1 

47T 



9r 



r%r 2 \ 

r 2 J 



cos 6 



(3* 



This density is then multiplied by the Jacobian and inserted into the Poisson solver with 
the proper boundary conditions. 

In the axisymmetric system (/ = r), the spherical radius satisfies r 2 = r 2 + z 2 . The 
physical domain consists of an outer sphere of radius r s2 is created in which an inner 
sphere of radius r a \ has been removed. The corresponding physical charge density is 
obtained via (122bh . 



1 

Air 



12 r% 



4- 



2r si r 



s2 



cos 9 



(39) 



The L2-norm of the error for both azimuthal and axisymmetric solvers is shown in 
Table [3l Here R = r%jT\ {oir S 2/r s i) is the ratio of the radius of the outer boundary to 
that of the inner boundary and N^ v is the number of uniformly spaced grid points in 
each direction. Figure H] displays the chosen MMS potential on both the physical and 
logical domains as obtained by the logical grid Poisson solver. Notice that the different 



The ACC-PIC Method 



16 




.2 D D.Z D.4 



(a) 



(b) 



Figure 4. Two-dimensional MMS potential plots as obtained computationally using 
and in physical (a) and logical (b) space. 



values of p^ inserted in the solver for the different coordinate systems should (and do) 
return the same potential on the grid. The results shown in Tables [2] and [3] show the 
expected second-order scaling with grid spacing, showing that our field solver is working 
correctly. 



6.3. Self-Fields and Momentum Conservation 

With a uniform grid, symmetrical particle shapes guarantee exact momentum 
conservation if the gridding is done correctly [3J . For a uniform rectangular grid, this is 
equivalent to having the self-electric forces exactly zero on each particle. (126bl) dictates 
that we use the logical grid electric field at the particle's position on the logical grid. 
The most obvious method to do this is by simply differencing $ with respect to £, 
(where we have defined a staggered mesh such that exists on vertices and $ and p^ 
exist on cell centers). One then interpolates E^ to the logical space particle position. 
Unfortunately, simple Id tests on this method of obtaining the logical electric field with 
a single particle reveal a non-zero self-field at the particle. 

We have therefore devised a second, less direct method of obtaining the logical 
electric fields at the particle position in which we solve for the electric field on the 
logical grid and interpolate this field to the particle position. As such, we calculate the 
electric field on the vertices, formulated on the logical grid as 



E 



and 



E v ■ 

h3 



2A£ 



2Ar) 



(40a) 



(40b) 



We then apply the proper extrapolation techniques such that the overall second-order 
accuracy of the system is upheld and calculate the physical electric fields on the vertices 
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using 



TV 
J i,3 



and 



'■■'I, 7r(./; l .,,/^-./;,, ; /^,). (41b) 

where we have transformed from the inverse Jacobian matrix, k^ u to the Jacobian matrix 
jfa,. Here, the superscript v on the components of the Jacobi matrix and its determinant 
signify that they are calculated on the vertices using a four-point average from their 
natural cell-centered locations. This technique for calculating the electric fields at the 
particles leads to exactly zero self-forces on a test particle in Id, but in 2d the Jacobi 
matrix components do not exactly cancel, leading to non-zero self-forces on the particles. 
Reference [17] provides a detailed explanation as to why this method leads to non-zero 
self-forces at the particle and a detailed comparison of this method with the direct 
interpolation method detailed above. We stress that, with a non- uniform grid, zero 
self-forces do not guarantee exact momentum conservation because of the presence of 
the inertial terms in (Tl5|) . which are often larger than the self- force terms |17j . 

6.3. 1. Logical Electric Fields at a Particle To test the self- forces on the particle in 2d, 
we have set up a system in which a single particle is at rest on a doubly-periodic grid 
with a neutralizing background such that we can interpolate both of the physical electric 
fields and required grid derivatives to the grid positions and retroactively multiply them 
to obtain the logical fields at the particle position. As mentioned above, we expect to 
have some small self-force in 2d by utilizing this method. Below we attempt to quantify 
these forces. 

In 2d, the logical electric fields are obtained from the physical fields using 
<9<3> dx u <9$ 



(42) 



such that 



Q£H Q£n Q x v 

Et = Jll E x + j21 Ey 

e v = jvi /•:•' • hi /-'"• 1 ] 

For a uniform grid, the electric fields at the particle are zero to machine precision. 
However, as soon as we allow one of the grid dimensions to become nonuniform, the 
fields at the particle in the nonuniform dimension are no longer zero, even for e g = 10~ 4 . 
We find that the fields at the particle position scale with second-order accuracy in grid 
spacing, and the magnitude of the field at the particle in the nonuniform dimension is 
dependent upon the magnitude of the grid nonuniformity parameter, e g . For example, 
using a grid with = N v = 32 where e g = 0.15 in x and holding y to be uniform 
(e g = 0), E*(£ p ,ri P ) « 10" 5 , whereas r] p ) is zero to machine precision. For e g = 0.1, 

same grid resolution. 
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(a) 



(b) 



Figure 5. Comparison of cold plasma oscillation field energies using a uniform grid and 
the nonuniform grids given by (|25a[) and (|25b[) for (a) a long run and (b) a zoomed 
section of the run showing an oscillation period of 2ir for all three grids. Here we 
have used quadratic (5*2) particle shape functions with N% = N v = 128, N ppc = 225, 
At = 0.025, and e g = 0.1 (for the nonuniform case). 



We have also tried the direct interpolation method, i.e. directly interpolating the 
logical electric field on the logical grid to the particle position. This method also 
scales with second-order accuracy in grid spacing, but the indirect method detailed 
above consistently provides smaller fields at the particle position. 

7. NUMERICAL RESULTS 

We have performed several benchmarking tests to validate the full 2D curvilinear 
coordinate PIC code. For all tests presented below, we use a nonuniformity parameter 
of e g = 0.1 (for both orthogonal and nonorthogonal grids). This leads to a ratio of 
areas of the largest cell to smallest of ~ 20 for the nonuniform orthogonal grid (I25aj) 
and ~ 4.4 for the nonorthogonal grid (I25bj) . The maximum skewness parameter S max 



is ~ 0.75 for the nonorthogonal grid with this nonuniformity parameter. We have done 
only high frequency tests and have accordingly assumed that the ions are an immobile 
background, providing charge neutrality in the equilibrium. We have normalized the 
equations so that the electron plasma frequency u pe is equal to unity. 

7.1. Cold Plasma Oscillations on a Square Physical Domain 

Figure |5] shows a comparison of the evolution of the electrostatic field energy, defined in 
2d as 

<E 2 >= l -j dx dy [{E x f + [Eyf] , (44) 
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for a cold plasma oscillation on the uniform, nonuniform but orthogonal, and 
nonorthogonal grids. For these tests, we take a cold distribution of electrons with a 
uniform density stationary ion background and initialize a perturbation in the system 
by perturbing the particle positions with respect to the uniform neutralizing background. 
We define the quantity e pcrt to be the size of this initial perturbation on the particle 
positions. The initial particle positions, whether assigned uniformly in the physical 
space or assigned according to the Jacobian in the logical space, lead to a perturbed 
charge density. This perturbation can dominate e per t if the latter is small enough. For 
these tests, the initial particle perturbation is directed at a 45° angle across the grid to 
check the effects of the interpolation in multiple dimensions on the data produced. In 
Figure we show both a long-time evolution and a zoomed section of the field energy 
from these same runs. We have used a perturbation of ep ert = (7.07e~ 5 , 7.07e~ 5 ), such 
that the magnitude of the perturbation is |ep e rt| = 1 x 10~ 4 , N% = N v = 128, the average 
number of particles per cell is N ppc = 225, and u pe At = 0.025. No measurable growth 
in the electrostatic field energy was observed during the course of these runs. We note 
that the period of the plasma oscillations observed in Figure |5] is very nearly 2ir (as is 
expected for u pe = 1) for all three grid choices. Further testing has revealed that the 
period of the plasma oscillations converges to the value of 2ir with second-order accuracy 
in At, as expected, and the electric field energy does not begin to decay even in very 
long runs. 



7.2. Cold Electron- Electron Two-Stream Instability 
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Figure 6. Comparison of cold electron-electron two-stream instability growth rates 
for uniform and non-uniform grids. Higher initial noise levels due in the nonuniform 
grid case provide a larger initial perturbation in the system, leading to the differences 
seen between the two curves. Here we have used quadratic particle shape functions 
with N% = N v = 128, N ppc = 225, At = 0.025 for all cases, with e g = 0.1 for the 
nonuniform orthogonal and nonorthogonal cases. 



Having chosen a physical grid such that x, y G [— 7r : 7r], the wave vector k is given 
by k = [1,1], such that, to generate a two stream instability at a 45° degree angle 
across the grid, we use the effective wave vector, k 45 ° = y/2 fa 1.414. In our normalized 
units, the maximum growth rate for the two stream instability of ~ 0.35 occurs at 
\k ■ vq\ m 0.63. As such, we have chosen an initial velocity parallel to k (in the physical 
space) with v= [0.314,0.314]. 

Figure M shows the time evolution of the electrostatic field energy for a uniform grid 
compared with the nonuniform grids given by f)25ap and (125b|) . The offset of the three 
curves is due to the larger initial perturbation due to the grid non-uniformity discussed 
above. Here we have used quadratic particle shape functions with = N v = 128, 
Np P c = 225, oo pe At = 0.025 for all cases, again with e g = 0.1 for the nonuniform 
orthogonal and nonorthogonal cases. The same cases have been performed using bilinear 
particle shape functions, revealing identical results. For all cases, the observed growth 
rate is ~ 0.35, which matches the theoretical prediction. 



7.3. Landau Damping 

In this section we study the effects of Landau damping with our 2D code, on a unit square 
physical grid. Since the our thermal distribution is taken to be Maxwellian, the particle 
velocities are isotropic in x, y, and as such we have set up a case in which the initial 
perturbation is only in x for simplicity. Figure [7] shows the time evolution of the Landau 
damping on our electrostatic field energy for the uniform, nonuniform orthogonal and 
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Figure 7. Comparison of Landau damping rates for a uniform grid with those obtained 
using the nonuniform, orthogonal and nonorthogonal square grids given by (12 5 a I) 
and (|25bj) . Here we have used N% = N v = 128, iV ppc = 400 per species, At = 0.025, 
Vth = 0.07, and quadratic particle shape functions for all cases, with e g = 0.1 for the 
nonuniform orthogonal and nonorthogonal cases. 



nonorthogonal square grids. Here we have used N% = N v = 128, iVpp C = 400, At = 0.025, 
and Vth = 0.07 for all cases. Note that the nonorthogonal grid more closely follows the 
damping rate of the uniform grid than does the nonuniform, orthogonal grid. The 
real oscillation frequency agrees very well with the theoretical value of 0.58 given by 
the Bohm-Gross dispersion relation [25] u 2 = u 2 e + 3k 2 v 2 h , where again we have taken 
u pe = 1 in accordance with our code normalizations. The average damping rate across 
the simulation time agrees to within 3% of the theoretical value of 0.124 as given by [25] 




for these parameters for the uniform and nonorthogonal grids, whereas the damping 
rate on nonuniform, orthogonal grid agrees to within 5%. We hypothesize that the 
lower damping rate observed for the nonuniform orthogonal grid case (I25bl) is due to 
the strong nonuniformity of the grid; the ratio of the largest to the smallest cell areas 
for a nonuniformity parameter e g = 0.1 is ~ 20. To check this hypothesis, we have also 
run cases in which we have used e g = 0.06 such that the largest to smallest cell area 
ratio is ~ 4.9. These cases give the same damping rates as the uniform grid case shown 
above. 

7.4- Cold Plasma Oscillations on a Concentric Annulus 





(45) 



As a final test of our entire method, we have set up a cold plasma oscillation on the 
concentric annulus grid of Figure (taking azimuthal symmetry). As an initial test, 
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Figure 8. Comparison of cold plasma oscillation field energies on an annular physical 
grid using an initial perturbation in r only and a combination of r and initial 
perturbation showing an oscillation period very close to 2n for both cases. Here we 
have used quadratic spline particle shape functions with = N 7] = 64, N ppc = 400, 
and At = 0.05. 

we perturbed the initial potential radially using 

<l = e p crtCOS ( 7r (^r))' ( 46 ) 

where we have used e pcrt = 10 -4 , T\ = 0.25, and T2 = 1.0. This initial perturbation 
satisfies homogeneous Neumann boundary conditions along the entire boundary of the 
system. 

In general, for Neumann boundary conditions, there is a subtlety, namely that 
the unit normal to the physical dommain does not map to the unit normal in the 
logical domain. Thus, Neumann boundary conditions involve the normal and tangential 
derivatives on the logical boundary. For this case, however, the Winslow coordinates are 
orthogonal, and this issue does not arise. See Sec. 6.2 for a discussion of the associated 
null space issues. 

With our input parameters, the ratio of the area of the largest grid cell to the 
smallest is ~ 16. This is a very simple test of the system, and is analogous to perturbing 
our system only in y on a rectangular grid. The time evolution of the electrostatic field 
energy for this test is shown in the blue curve of Figure M for a 64 x 64 physical grid 
with u pe At = 0.1, iVpp C = 225, and quadratic spline shape functions. The period of 
oscillation of the field energy for this set of initial conditions is observed to be ~ 6.3 for 
this case, meaning that the frequency is indeed unity with an error that scales as A£ 2 , 
as per our code normalizations. 
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As a more challenging case, we then perturbed the initial potential using 

$ = e pcrt cos ^tt _^ J ) cos 9, (47) 

where we have again used e pert = 10~ 4 , r\ = 0.25, and r 2 = 1.0. The red curve in 
Figure [5] shows the time evolution of the electrostatic field energy for a 64 x 64 physical 
grid with oj pe At = 0.05, N ppc = 400, and quadratic spline shape functions. We have used 
more resolution in this particular case in order to resolve more fully the complicated 
features of our initial potential. Again, the measured period of oscillation is ~ 6.3 (the 
difference between the two curves is approximately 0.085%); giving frequency unity, 
with errors scaling as A£ 2 . 




X 



(c) uj pe t = 3.6 (d) oj pe t = 4.8 

Figure 9. Snapshots of one period of the evolution the potential of a cold plasma 
oscillation on the annular physical grid. Here we have perturbed the inital potential 
using (|47|) . 

Figure |9] shows the time evolution of one period of the potential on the physical 
grid using the initial perturbation as given by ( I47jl for the cold electrostatic plasma 
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oscillation test described above. Figures M^&) and (b) are from one half of the plasma 
period and Figures M^c) and (d) are from the next. Notice the exact reversal of the 
potentials between the two halves of the plasma period. Reflecting boundary conditions 
for particles at the non-periodic parts of the boundary were not required for this case, 
because the particle excursions for a cold plasma oscillation are so small. 

8. CONCLUSIONS 

In this study, we have set out to develop a new approach to the PIC method, in which the 
key components, the mover, the field solve, charge accumulation, and field interpolation, 
are carried out on a uniform grid in the logical domain, the unit square, and mapped 
to an arbitrary boundary conforming grid on an arbitrary physical domain. We have 
focused on a 2d, electrostatic model as a proof of principle. 

To demonstrate our method, we used analytically prescribed grids and grids 
obtained by Winslow's grid generation technique [15] to generate a boundary-fitted grid, 
in the latter case by solving a set of coupled elliptic equations using a nonlinear Newton- 
Krylov solver. The generated grids are then mapped onto the logical grid through the 
use of the mapping x(£ ) and its inverse, such that the PIC components can be 

run on the logical grid using these mappings. 

We have derived the logical grid macroparticle equations of motion based on a 
canonical transformation of Hamilton's equations from the physical domain to the 
logical. The resulting nonseparable system of equations using an extended form of the 
semi- implicit modified leapfrog integrator [21], which we have shown to be symplectic 
for a system of arbitrary dimension and second-order accurate in time if symmetrization 
leading to a time-centered discretization of the macroparticle update equations is 
utilized. If the field solve is performed just after the first step of this integrator, it 
needs to be done only once per time step. 

In order to obtain the electrostatic fields on the logical grid, we have constructed a 
generalized curvilinear coordinate formulation of Poisson's equation which is discretized 
conservatively on the logical grid using a staggered mesh. Field boundary conditions 
are applied in such a way as to produce a symmetric operator matrix which we then 
solve using a conjugate gradient solver. 

Our formulation of the curvilinear coordinate Poisson equation requires the logical 
grid charge density, which allows us to accumulate the charge from the particles directly 
onto the uniform, square logical grid using standard particle shape functions rather 
than the more complicated, weighted shape functions which must be used if the charge 
is accumulated on a nonuniform physical grid. Furthermore, the particle equations of 
motion require that the derivative of the electrostatic potential on the logical grid be 
obtained at the particle positions for the update of the particle momentum. These 
logical electric fields are interpolated to the particle positions on the logical grid using 
the symmetric particle shape functions which have been slightly modified from the 
standard shape functions used in the charge accumulation process in order to account 
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for our choice of a staggered mesh. All validation tests performed have shown that 
the code produces correctly the physics for complicated interior meshing, as well as for 
complicated domain shapes. 

Albeit at a proof-of-principle level, our code has shown that an arbitrary coordinate 
PIC method is a an accurate and efficient alternative to current methods for simulating 
plasma systems with complex domains. Future work will focus on extending the 
method to 3d and extending to electromagnetic systems. Further, work is required 
on parallelization techniques to handle more efficiently the particle push and charge 
accumulation stages of our method. Finally, our method can be coupled to a moving 
mesh algorithm to more accurately follow the dynamic evolution of the plasma system. 
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Appendix A. Differential Geometry Notation 

This Appendix presents the reader with a coherent overview of the relations between 
Cartesian and curvilinear coordinates which form the basis of this paper. 

Let the values x a , a = 1, • ■ ■ , n be the Cartesian coordinates of the vector x. The 
coordinate transformation x(£ ) defines a set of curvilinear coordinates £ a , ■ ■ ■ , £ n in the 
domain X n . We define the Jacobi matrix of this transformation 



Conversely, we can also think of this transformation as a mapping of £ to x, £,(x). 
Defining the inverse of the matrix j a/ 3 as 




n, 



(A.l) 



and its Jacobian 



J(i) = det (j) . 



(A.2) 




a, (3 = 1, ■ ■ • 



(A.3) 



we can write its Jacobian 



K{x) = det (fc) 



1 



(A.4) 



J 
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Simple linear algebra tells us that 

f)r a r)fP 

^ = W^ = *» (A - 5) 

therefore jk = I. 

Given a Euclidean metric on the physical space, we have 

dx^dx^ = ——dCde = gafidFdZP, (A.6) 
and the covariant metric tensor is defined as 

9ap{i) = d £ ad £f) > «,/3,7 = l,---,n. (A.7) 

From flA.ip . we see g a p = j 7a j 7 p, which is simply gr cov = j T j. Likewise, the contravariant 
metric tensor is based upon the inverse Jacobian matrix, k: 

9 = = !,■••,«, (A.8) 

thus g con = kk T . 

Finally, by multiplying the covariant and contravariant metric tensors and applying 
the identity jk = I, we can prove that the covariant and contravariant metric tensors 
are in fact inverses of each other: 

<7cov<7 C ° n = fjkk T = j T k T = (kj) T — I T — I. (A.9) 

In two dimensions, we can easily convert from the covariant to the contravariant metric 
tensor using the following equation: 

g ap = 93-«,3-^ a ,P = l,2, (A.10) 

fi'cov 

where g cov = J 2 is the determinant of the covariant metric tensor, and g cov = where 
gcon _ j£i j g determinant of the contravariant metric tensor. Likewise, we can shift 
from the contravariant to the covariant metric tensor using 

gap = (-i) a ^g cov g 3 - a > 3 ^, a,/? = 1,2. (All) 

Appendix B. Derivation of Curvilinear Coordinate Poisson Equation 

We begin by rewriting ([3]) as the divergence of the gradient of the potential in generalized 
coordinates: 

1 1 d <9$ 

where / is a geometry factor allowing us to switch between azimuthal and axially- 
symmetric systems. For instance, if we want the physical coordinate system to be x, y 
with = 0, we set / = 1. For r, z coordinates with = 0, we set f — r. 
In matrix notation we can write V<3? as 

™ = 

= A m Vr (B.2) 
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where A m are the covariant components of the 2D vector 

A = A m Vr = A 1 V£ + A 2 Vr ] . (B.3) 
Likewise, we can represent A by its contravariant components, A m : 

A = A l Vr) x z + A 2 z x Vf, (B.4) 
such that we can write 

AiVf + A 2 Vr] = A x Vr] x z + A 2 zx Vf . (B.5) 

We can now find the direct relationship between the covariant and contravariant 
components of A. For example, taking the inner product of flB.51) with Vrj x z we 
find: 



A 1 (V77 x 5) • V£ = A|V£| 2 + A 2 Vr/ • Vf . 

By (TO, 

• vr = <A 

allowing us to rewrite flB.61) as 

A^ 11 + A 2 g 12 = A 1 (V77 xf)-V( = A X V£ • V?? x z. 

We note that 

<9£ dr/ dt; dr/ 
dx dy dy dx 
= det(fe) 
1 

~ J' 

so we can write the left-hand side of (1B.8[) as 

A 1 



V£ • V?7 x z 



12 



Similarly, 



A 1 g 11 + A 2 g 
A l9 21 + A 2 g 22 = ^. 



Returning now to ( IB. If) and using ( 1B.2I) and ( 1B.4j) and the identity 
V • (cv) = Vc ■ v + cV ■ u, 
we can write the Laplacian operator as 

ly • /V$ = V • A 
/ 

= i [VifA 1 ) ■ (Vrj xz) + V(fA 2 ) ■ (z x VO] . 
Using V£ • (V£ x z) — 0, (IB. 12j) can be written as 



(B.6) 
(B.7) 
(B.8) 



(B.9) 

(B.lOa) 

(B.lOb) 
(B.H) 



- V • /V$ 



^Vtv^ + I^v^xvo' 



/ ^ 



. (B.12) 
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By (IB.9I) . we can rewrite (IB. 12ft as 

Now converting the contravariant components of A to their covariant formulations using 
( IB. 21) and ( IB. 81) . we can write the final form of the Poisson equation in curvilinear 
coordinates: 

= (B.14) 



where p = p x is the physical charge density. 
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